Compute stress and strain#
In this demo we show how to compute stress and strain
import dolfin
try:
from dolfin_adjoint import Constant, DirichletBC, Function, project
except ImportError:
from dolfin import DirichletBC, Constant, project, Function
import pulse
from fenics_plotly import plot
gamma_space = "R_0"
geometry = pulse.Geometry.from_file(pulse.mesh_paths["prolate_ellipsoid"])
# geometry = pulse.geometries.prolate_ellipsoid_geometry(mesh_size_factor=3.0)
geometry.mesh.coordinates()[:] *= 0.1
# breakpoint()
activation = Function(dolfin.FunctionSpace(geometry.mesh, "R", 0))
2024-03-27 07:47:47,073 [499] INFO pulse.geometry_utils:
Load mesh from h5
matparams = pulse.HolzapfelOgden.default_parameters()
material = pulse.HolzapfelOgden(
active_model=pulse.ActiveModels.active_stress,
T_ref=1.0, # Total active stress should be activation * T_ref
eta=0.2, # Fraction of transverse stress
activation=activation,
parameters=matparams,
f0=geometry.f0,
s0=geometry.s0,
n0=geometry.n0,
)
# LV Pressure
lvp = Constant(0.0)
lv_marker = geometry.markers["ENDO"][0]
lv_pressure = pulse.NeumannBC(traction=lvp, marker=lv_marker, name="lv")
neumann_bc = [lv_pressure]
# Add spring term at the base with stiffness 1.0 kPa/cm^2
base_spring = 1.0
robin_bc = [
pulse.RobinBC(value=Constant(base_spring), marker=geometry.markers["BASE"][0]),
]
def fix_basal_plane(W):
V = W if W.sub(0).num_sub_spaces() == 0 else W.sub(0)
bc = DirichletBC(
V.sub(0),
Constant(0.0),
geometry.ffun,
geometry.markers["BASE"][0],
)
return bc
dirichlet_bc = [fix_basal_plane]
# Collect boundary conditions
bcs = pulse.BoundaryConditions(
dirichlet=dirichlet_bc,
neumann=neumann_bc,
robin=robin_bc,
)
# Create the problem
problem = pulse.MechanicsProblem(geometry, material, bcs)
# Solve the problem
pulse.iterate.iterate(problem, lvp, 15.0, max_nr_crash=100, max_iters=100)
pulse.iterate.iterate(problem, activation, 60.0)
2024-03-27 07:47:47,144 [499] INFO pulse.iterate: Iterating to target control (f_35)...
2024-03-27 07:47:47,145 [499] INFO pulse.iterate: Current control: f_35 = 0.000
2024-03-27 07:47:47,145 [499] INFO pulse.iterate: Target: 15.000
2024-03-27 07:47:47,269 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,270 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,271 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,271 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,272 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,272 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,273 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,273 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,274 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,275 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,275 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,276 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:47,276 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:47:48,262 [499] DEBUG UFL_LEGACY: Blocks of each mode:
99 full
2024-03-27 07:47:48,337 [499] DEBUG UFL_LEGACY: Blocks of each mode:
2024-03-27 07:47:48,544 [499] DEBUG UFL_LEGACY: Blocks of each mode:
18 full
2024-03-27 07:48:51,703 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,704 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,705 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,705 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,706 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,706 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,707 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,707 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,708 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,708 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,709 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,709 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,710 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:48:51,893 [499] DEBUG UFL_LEGACY: Blocks of each mode:
10 full
2024-03-27 07:48:51,968 [499] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
2024-03-27 07:48:52,162 [499] DEBUG UFL_LEGACY: Blocks of each mode:
3 full
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_DTOL.
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_FNORM_NAN.
*** Warning: PETSc SNES solver diverged in 1 iterations with divergence reason DIVERGED_FNORM_NAN.
2024-03-27 07:49:01,713 [499] INFO pulse.iterate: Iterating to target control (f_22)...
2024-03-27 07:49:01,714 [499] INFO pulse.iterate: Current control: f_22 = 0.000
2024-03-27 07:49:01,714 [499] INFO pulse.iterate: Target: 60.000
*** Warning: PETSc SNES solver diverged in 0 iterations with divergence reason DIVERGED_FNORM_NAN.
([Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1202),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1266),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1391),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1454),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1533),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), MixedElement(VectorElement(FiniteElement('Lagrange', tetrahedron, 2), dim=3), FiniteElement('Lagrange', tetrahedron, 1))), 1597)],
[Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1200),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1264),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1389),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1452),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1531),
Coefficient(FunctionSpace(Mesh(VectorElement(FiniteElement('Lagrange', tetrahedron, 1), dim=3), 1), FiniteElement('Real', tetrahedron, 0)), 1595)])
V = dolfin.VectorFunctionSpace(geometry.mesh, "CG", 1)
u, p = problem.state.split(deepcopy=True)
# u_int = interpolate(u, V)
# mesh = Mesh(geometry.mesh)
# dolfin.ALE.move(mesh, u_int)
# Get the solution
# u, p = problem.state.split(deepcopy=True)
# dolfin.File("u.pvd") << u
W = dolfin.FunctionSpace(geometry.mesh, "CG", 1)
F = pulse.kinematics.DeformationGradient(u)
E = pulse.kinematics.GreenLagrangeStrain(F)
# Green strain normal to fiber direction
Ef = project(dolfin.inner(E * geometry.f0, geometry.f0), W)
# Save to file for visualization in paraview
# dolfin.File("Ef.pvd") << Ef
plot(Ef)
2024-03-27 07:49:03,947 [499] DEBUG UFL_LEGACY: Blocks of each mode:
2024-03-27 07:49:04,351 [499] DEBUG UFL_LEGACY: Blocks of each mode:
1 full
<fenics_plotly.fenics_plotly.FEniCSPlotFig at 0x7f21b70c46a0>
P = material.FirstPiolaStress(F, p)
# First piola stress normal to fiber direction
Pf = project(dolfin.inner(P * geometry.f0, geometry.f0), W)
# Save to file for visualization in paraview
# dolfin.File("Pf.pvd") << Pf
plot(Pf, colorscale="viridis")
2024-03-27 07:49:05,198 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,198 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,199 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,199 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,200 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,201 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,201 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,202 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,202 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,203 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,203 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:05,446 [499] DEBUG UFL_LEGACY: Blocks of each mode:
1 full
<fenics_plotly.fenics_plotly.FEniCSPlotFig at 0x7f21b710a470>
T = material.CauchyStress(F, p)
f = F * geometry.f0
# Cauchy fiber stress
Tf = dolfin.project(dolfin.inner(T * f, f), W)
# Save to file for visualization in paraview
# dolfin.File("Tf.pvd") << Tf
plot(Tf, colorscale="jet")
2024-03-27 07:49:06,358 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,358 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,359 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,360 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,360 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,361 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,361 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,362 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,362 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,363 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,363 [499] INFO UFL_LEGACY: Adjusting missing element cell to tetrahedron.
2024-03-27 07:49:06,613 [499] DEBUG UFL_LEGACY: Blocks of each mode:
1 full
<fenics_plotly.fenics_plotly.FEniCSPlotFig at 0x7f21b7e47a90>